JYFL-5/00 
hep-ph/00103T9| 
20 October, 2000 



TRANSVERSE ENERGY FROM MINUETS 
IN ULTRARELATIVISTIC NUCLEAR COLLISIONS: 
A NEXT-TO-LEADING ORDER ANALYSIS 



K.J. Eskola|J' fe and K. Tuonanier^' 



a Department of Physics, University of Jyvdskyld, 
P.O.Box 35, FIN-40351, Jyvdskyld, Finland 

b Helsinki Institute of Physics, 
P.O.Box 9, FIN-00014, University of Helsinki, Finland 



Abstract 

We compute in next-to-leading order (NLO) perturbative QCD the amount of trans- 
verse energy produced into a rapidity region AY of a nuclear collision from partons 
created in the few-GeV subcollisions. The NLO formulation assumes collinear fac- 
torization and is based on the subtraction method. We first study the results as a 
function of the minimum transverse momentum scale and define and determine the 
associated .fT-factors. The dependence of the NLO results on the scale choice and on 
the size of AY is also studied. The calculations are performed for GRV94 and CTEQ5 
sets of parton distributions. Also the effect of nuclear shadowing to the NLO results 
is investigated. The main results are that the NLO results are stable relative to the 
leading-order (LO) ones even in the few-GeV domain and that there exists a new kine- 
matical region not present in the LO, which is responsible for the magnitude of the 
.fT-factors. The dependence on the size of Ay is found to be practically linear within 
two central units of rapidity, and the scale dependence present at RHIC energies is 
seen to vanish at the LHC energies. Finally, the effect of NLO corrections to the de- 
termination of transverse energies and multiplicities in nuclear collisions is discussed 
within the framework of parton saturation. 
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1 Introduction 



The heavy ion community experiences currently a new exciting period as the first 
data in ultrarelativistic heavy ion collisions at BNL-RHIC are being analysed and the 
first physics results are coming out. These include the observation by the PHOBOS 
collaboration of an increased activity in Au+Au collisions at \fs = 56 and 130 AGeV 
[|1]: the number of charged particles per participant pair clearly increases from that 
in Pb+Pb collisions at the CERN-SPS at y/s = 17.3 AGeV, and from pp at y/s = 
50 . . . 200 GeV. This lends support to a description of initial particle production in high 
energy heavy ion collisions in terms of binary collision mechanisms. The first results 
from the STAR collaboration have in turn shown an increase of elliptic flow in Au+Au 
at \/s = 130 AGeV, observed as an enhanced asymmetry in the azimuthal distributions 
of particles f2j. This points towards collectivity, pressure and thermalization in the 
produced system. 

A complete description of a collision of two heavy ions at high energy is clearly an 
extremely challenging task: elementary QCD-quanta, gluons, quarks and antiquarks, 
are expected to be produced abundantly enough, so that they thermalize locally and 
form a new, deconfined phase of matter, the Quark Gluon Plasma (QGP). After its 
production the system will experience an expansion stage during which it will cool 
down and go through a phase transition from the QGP to a hadron gas. Eventually, 
decoupling sets in and the hadrons fly independently to the detector, possibly decaying 
into other particles before that. 

Due to our ignorance of how to solve nonperturbative QCD, calculations to de- 
scribe such a complicated spacetime evolution are very difficult to do from truly first 
principles. Instead, the evolution can be effectively modelled in terms of relativistic 
hydrodynamics || [|, which relies on the assumption of a locally thermal system. The 
benefit in such an approach, in addition to simplicity, is a proper treatment of the phase 
transition and inclusion of collective effects (pressure and flow). The initial conditions, 
however, have to be given from outside. Alternatively, one may try to describe the 
strongly interacting system first as a cascade of partons H |6| and then as a hadronic 
cascade [0, ||, or, as recently suggested, as a combination of the hydrodynamical ap- 
proach and a hadron cascade ||. 

Common to the approaches mentioned above, the initial conditions, i.e. the energy 
densities and number densities of initially produced quarks and gluons, along with 
the formation times, need to be known as precisely as possible. As also the degree of 
thermalization of the QGP |TI], |IT], |T^] depends heavily on the initial conditions, it is 
vitally important to study them in detail. In this paper we will do this by formulating 
and computing the initial transverse energy production from minijets in central AA 
collisions to next-to- leading order (NLO) in perturbative QCD (pQCD). 

Initial particle production in AA collisions at collider energies has typically been 
described with two-component models where pQCD accounts for the hard and semihard 
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scatterings of partons and a phenomenological model for the soft particle production 
l3j 0, [15|, 16|. At sufficiently high cms-energies, partons produced with transverse 
momenta in the semihard domain, px > po ~ 1 ... 2 GeV, become the dominant 
source of energy deposited in central rapidities [|13|, [L7|, [18j . Within the leading twist 
framework the inclusive cross sections for production of semihard partons, along with 
the initial E? carried by them, are computable in terms of collinear factorization |19 
of pQCD 0, Individual binary scatterings of partons (or nucleons) are assumed 
to be independent of each other and nuclear effects can be included by using nuclear 
parton distributions [pOfl . Computation of the pQCD component thus leans on the 
fact that production of quarks and gluons is reliably described by pQCD at large 
transverse momentum scales, where it is possible to actually observe individual jets in 
pp collisions. The main uncertainties are related to the determination of the semihard 
scale po defining the perturbative component, and to the higher order contributions 
to the parton cross sections. To reduce the latter source of uncertainty is the main 
motivation for the present NLO study. 

The studies |13| , |T4t |T5| , 17| are in practise extensions from the pp {pp) physics in the 
sense that the smallest transverse momentum included in the perturbative component, 
Po, is constrained by data from hadronic physics and kept constant in the extrapolation 
to AA. Recently, however, ideas of parton saturation J[R |21| have been revived for 
an effective estimation of the initial gluon and quark production based on the pQCD 
component [ Z2"|, alone. The perturbative production of partons is extendend down 
to a minimum scale p = p sat , defined as the scale at which the produced partons, each 
of which is of the size vt/pq, effectively fill the available transverse area ttR\. This is the 
simplest way of defining a saturation criterion in central AA collisions. For A ~ 200, 
this procedure leads to p sat ~ 1 GeV for RHIC energies -y/i = 56 . . . 200 AGeV and 
p sat ~ 2 GeV for the LHC energy ^Js = 5500 AGeV [2^]. The criterion which deter- 
mines when the system becomes saturated, may contain nonperturbative, even phe- 
nomenological features. Terminating the pQCD computation at p Q = p sa ,t(y/s,A) thus 
gives an estimate of the total initial parton production. The effective formation time of 
the minijet (gluon) plasma system is then obtained as r ~ l/p s &t- The determination 
of p sat in the pQCD+saturation approach is a dynamical procedure which is affected 
by the NLO contributions to partonic cross sections. Again, reducing the uncertainty 
related to the higher order contributions in minijet production is the general goal in 
this paper. 

Another promising approach for describing initial gluon production is that of clas- 
sical fields 
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The gluon fields originate from colour sources moving along the light 
cones and, assuming longitudinal boost-invariance, an effective dimensionally reduced 
Hamiltonian can be formed and the gluon production computed through a lattice sim- 
ulation [25|| . Interestingly, also in these classical field computations gluons with fairly 
hard momenta E T /N ~ 3Q S at, where Q sat is a similar saturation scale as p sat above, 

As this approach should coincide with 



are important for the initial production |26|, |27 



2 



the perturbative one in the limit of high transverse momentum of produced gluons, 
it is important to reduce the uncertainties related to the possible higher-order contri- 
butions in the pQCD component. This gives further motivation to perform the NLO 
computation of initial Et we present in this paper. 

In our study the NLO evaluation of the produced initial transverse energy in central 
(6 = 0) AA collisions is carried out by computing the first E T moment <tqcd(Et)ay, P0 
of the perturbative Et distribution of minijets in a single proton-proton collision within 
the framework of collinear factorization. Nuclear collision geometry is taken into ac- 
count through the standard nuclear overlap functions Taa{^>) and the initial Et is 
obtained as E T = Taa&qcd{Et) [0- The parameter p defines the extent of the per- 
turbative treatment: we include all hard scatterings of partons in which at least an 
amount of 2p of transverse momentum is created. Note that p ^> Aq CD is kept as 
an external parameter, allowing for different phenomenological ways to determine its 
value. 

We emphasize that although the initial Et in a rapidity acceptance window Ay is 
a physical quantity, which thus can be computed in an infrared safe manner, it is not 
a direct observable. This is due to the interactions of initially produced gluons and 
quarks which are expected to thermalize the QGP system rapidly. Due to the pressure 
P formed, the system does work (PdV) during the expansion stage, and some of the 
initially stopped energy is lost. In ref. |22) it was estimated that for central collisions 
of A ~ 200 the measurable final state E T is £^ itial /3 at the highest RHIC energy and 
£^ nitial /5 at the LHC energy. 

In practise, the NLO computation is performed by appealing to a subtraction algo- 
rithm |3(J], applicable for a wide class of problems requiring a NLO evaluation of 
a quantity fulfilling the requirements of infrared safeness. The required 0{a z s ) matrix 
elements have been evaluated in |]3~If . The current paper is a sequel to our previous 



work [Q, where <Jqcd(Et)ay, P0 was computed numerically and the results studied as 
a function of the parameter po and the factorization/renormalization scale choice. The 



main result of |32[, that the deviation of the NLO results from the LO ones are not 
dramatically increasing as po is taken down from 10 GeV to 1 ... 2 GeV, remains also 
here. In |32| it was noted that a new and important kinematical perturbative region 
has to be included in the NLO computation. Here we study the contributions from this 
region in more detail. The relation to another recent NLO computation of the initial 



Et, Ref. |33fl , along with the differences from our formulation, were also discussed in 



52]. Now we extend the NLO study to explore the sensitivity of the results on the 
choice of parton distributions. In addition to the dependence of the results on p and 
on the scale choice, their dependence on AY is explicitly shown and the important 
implications discussed. Also the effects of nuclear parton distributions |2(| are stud- 



ied. Finally, as an application of the NLO computation, charged particle production 
for the nucleus-nucleus collisions currently under way at RHIC as well as for those 
to take place in the future at the LHC/ALICE are briefly considered by applying the 
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pQCD+saturation model 



2 Theoretical basis of the NLO formulation 

In the framework of independent parton-parton scatterings, the average Ex produced 
initially into the rapidity region Ay in an AA collision at an impact parameter b can 
be computed as fIB] 

E^ A (h, v^, Po ) = T Aj4 (b) ( T QCD ( J E T ) AyiP0 , (1) 



where T AA (h) is the standard nuclear overlap function ||13|| . We will focus on cen- 
tral collisions of large nuclei only, for which, using the Woods-Saxon nuclear profiles, 
Taa(0) « A 2 /ttR a . Numerically, T AA (0) « 30/mb for A ~ 200. In what follows, we 
will focus on the formulation and computation of <Jqct>(Et)ay, P0 , the first moment of 
the minijet Et distribution in pp collisions. Regarding the Et distribution in an AA- 
collision, it is the first and also the second moments of the pp-\eve\ distributions that 
play the key roles rather than the distributions themselves [l^ . 



Calculation of physical cross sections in NLO pQCD is a nontrivial task even if the 
scattering amplitudes are known at the desired order. For infrared safe quantities, 
however, a systematic method allowing for the cancellation of the soft and collinear 



singularities in the loop and bremsstrahlung contributions exists [BG, 34]. In this section 



we will briefly review this algorithm and thoroughly describe its implementation. 
2.1 Subtraction method 

To order we need to take into account processes where two incoming hadrons pro- 



duce either two or three final state partons. Keeping the notations of [50], the incoming 
partons are labelled by A and B and the outgoing partons as 1, 2 and 3. For the two- 
parton final state the appropriate kinematical variables are yi,y2,PT2 and 02- Trans- 
verse momentum conservation determines pri = Pt2 = Pt and <pi = + as we do not 
include any intrinsic transverse momentum. The momentum fractions xa and Xb of the 
incoming partons are determined by the conservation of energy and longitudinal mo- 
mentum, resulting in = ^=(e ±?a + e ±m ). Similarly, for the three-parton final state 
p T i = — (pt2 + Pts) and the suitable kinematical variables are yi, y 2 , V3,Pt2,Pt3, 4>2 
and 63. The fractional momenta become now i 12 = En. e ±yi _|_ Pzg e ±y2 _j_ Ex2. e ±y-A_ A n 

' V s V s V s 

inclusive hard cross section can be written in a general form 

d[PSh d[PSh S2{pi>1,p " ) + J d [ p5 ] 3 d^^'^) 

= J[ 2 ^2]+/[2^3], (2) 
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where it is implicit that the integrations take place in 4 — 2e spacetime dimensions. 
To condense the notation, we have denoted the 2 — > 2 and 2^3 differential partonic 
cross sections as 

da 2 ^ 2 da 2 ^ 2 



d[PS} 2 dp T dy x dy 2 t 
da 2 ^ 3 da 2 ^ 3 



(3) 



(4) 



d[PS} 3 dp T 2dpT3dy 1 dy 2 dy : 

which also defines our notation for the phase space volume elements d[PS]2 and d[PS]z. 

The measurement functions S 2 (pi,p 2 ) and S 3 (pi,p 2 ,p 3 ) which depend on the four- 
momenta of the final state partons, define the physical quantity the cross section of 
which is to be computed. Such a quantity can be e.g. a differential one-jet cross 
section , , L |2~R 1311, or a two-iet cross section -j — , da , , I p [BHI for production of 

apTdy'H l| I' I |l ' •> dpxi dpT 2 dyi a«/2 1 ti LCi_|J 1 

observable jets with a jet cone radius R. In our case, S 2 and S3 will be designed to give 
ZE^hypo' ^ e differential E T distribution of minijets which fall into a given rapidity 
acceptance window AY and which originate from perturbative collisions where at least 
an amount 2p of transverse momentum is produced. 



Let us briefly recapitulate the procedure presented in detail in |30| for the cancellation 
of singular contributions in the above cross section. The 2^2 contribution can be 
written under collinear factorization |T^] as 

x(\M(a A + a B -> a 1 + a 2 )\ 2 }S 2 (p^pi 2 l ). (5) 

The incoming parton flavors are denoted by a a and a B whereas the outgoing parton 
distributions are denoted by a\ and a 2 . The matrix element squared and summed 
over final spin and color and averaged over initial spins and colors is evaluated in 
4 — 2e dimensions and contains the ultraviolet renormalization via the MS scheme. The 
distributions f(a,x) are the modified parton distributions defined in the MS scheme. 
Similarly for the 2^3 part one has 

d(j2 ^ 1 at / W / ^ PT2PT3 1 , f , 

x—fB{aB,XB)(\M{a A + a B ^ a x + a 2 + a-^S^p^p^p 1 ^), (6) 
x B 

where parton 3 has been identified as the one having smallest transverse momentum, 
canceling a factor 3 in the original statistical prefactor 1/3!. 

The two terms, 7 [2 — > 2] and I [2 — * 3] both contain divergent contributions pro- 
portional to 1/e and 1/e 2 . All these divergent contributions cancel remarkably against 
each other provided that the functions S 2 and S3 in the equation above, the so called 
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measurement functions, are infrared safe. This requirement can be stated more pre- 
cisely: First, as two of the outgoing partons become collinear, S3 should reduce to S 2 
i.e. 

S 8 (p?,(l-A)j&Ap£) = SiCpM, 

s 3 ((i-x)pl,p^Xpl) = sMrf), 

5 3 (AK,(l-AK,^) = ftCpf.rf) (7) 

for < A < 1. Second, as any of the partons becomes collinear with one of the beam 
momenta p A , Pb ft is required similarly, that 

$i(Pi,P2,Ap£) = S 3 (jft,j$, \p B ) = S 2 (pt,p%), 
S 3 (PiApa,P2) = Ss(Pi, Ap£,p£) = S 2 (pi,P2), 

Ss(k%,rt,rt) = s 3 (\p B ,pip») = sM,p 2 l ) (8) 

Assuming the infrared safety of the measurement functions, the calculation based 
on the subtraction method then proceeds by first separating the singular factors 
in 2 — > 3 matrix element from one another. As an initial step, the matrix element is 
decomposed as 

(|M| 2 > 2 _3 = (\M\ 2 ) A + (\M\ 2 ) B + (\M\ 2 ) 1 + (\M\ 2 ) 2 (9) 

in such a way that each term (\A4\ 2 ) n ~ 1/(P3 • Pn), so ft contains a soft singularity 
for parton 3 along with the singularity of parton 3 becoming collinear with parton n, 
where n = A,B,1,2 with A, B referring to the incoming partons and 1,2 to the two 
other outgoing partons. The imposed constraint ptz < Pti,Pt2 guarantees that other 
singularities do not occur. Now the 2 — ^ 3 part of the cross section can be written as 

J[2 -»• 3] = J[2 -> 3] A + I[2 -> 3] s + I[2 -» 3]i + J[2 -> 3] 2 . (10) 

Each of the four terms featured in the above equation can be treated independently 
and in a similar manner. For example, the last term becomes 

jp -* 3], = f§^4psu = fdipsu fy^f^^\ (id 

J d[PS\ 3 J p T3 [cosh(|/2 - ?/ 3 ) - cos(0 2 - 03j 

where the divergent factor l/(p 2 ' P3) h as been expressed in terms of the integration 
variables as p 2 ■ P3 = pT2PT3[cosh(y 2 — 2/3) — cos(0 2 — 03)], and F 2 is a complicated 
function containing the parton distributions, the measurement function S3 and certain 
parts of the (4 — 2e dimensional) squared matrix elements due to the decomposition in 
Eq. ©. 

The aim is now to decompose da^ 3 /d[PS]s into terms that are divergent but simple 
in a way that the integrations over the phase space of the third parton can be performed 
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analytically and terms that are finite as e — > but not so simple. This goal is achieved 
by inserting zero in the following clever manner: 

F 2 (yi,PT2,y2,4>2,PT3,y3,4>3) = 

F2(yi,PT2,y2,<fa,Pra,V3,<l>a) ~ F 2 (yi,p T 2i 2/2, 02, 0, y 3 , fa)9(pT3 < PT2/2) 

-F 2 (yi,PT2, 2/2, 4>2,PT3, 2/2, (j> 2 ) + F 2 (yi,PT2, 2/2, <f>2, 0, y 2 , $2)6 {PT3 < PT2/2) 

+F 2 (yi,PT2,y2,<fr2,0,y3,fo)0(pT3 < PT2/2) 
+F 2 (yi,PT2, 2/2, 02, PT3, 2/2,</»2) ~ ^2(2/1, PT2, 2/2, 02, 0, 2/2, 02)0(PT3 < ?W2). (12) 

The first four terms on the right-hand side combine to make up the term J [finite] 2 
which is perfectly finite as e — > as both the soft singularity px3 = and the collinear 
singularity 3 || 2 have been subtracted. The fifth term, containing the soft singularity, 
generates a term /[soft] 2, and the remaining two terms then form the term / [collinear] 2 
which, as the name suggests, contains the collinear singularity only. Thus 

I [2 -> 3] 2 = /[finite] 2 + /[soft] 2 + / [collinear] 2 . (13) 

Due to the decomposition (H) none of the above terms contain any other singularities, 
but these are then confined to other three terms I [2 — > 3]i, I [2 — > 3} a and I [2 — > 3]b 
defined and treated similarly. The procedure outlined here has therefore at this stage 
produced us the decomposition 

^[2^3] = Yl J [2^3] n 

n=A,B,l,2 

= J2 U[finite] n + /[soft] n + /[collinear] n }. (14) 

n=A,B,l,2 

The terms labelled 'finite' above, are very complicated but finite and can be evaluated 
numerically. In the singular terms /[soft] 2 and / [collinear] 2 the integrations over p^ 3 , 
03 and y 3 can be performed analytically due to the simple structure of those terms: the 
function F 2 is evaluated in each of them either on a soft or collinear limit. After these 
integrations are carried out, the resulting terms have 2^2 kinematics and contain 
several 1/e and 1/e 2 parts. Some of these cancel against identical (apart from the 
sign) terms in I [2 — > 2], and the remaining parts cancel among themselves beautifully, 
leaving a perfectly finite contribution behind - provided that the measurement functions 
S3 reduce to S 2 as required by Eqs. (0) and (|]). The somewhat complicated structure 
of the cancellation of the 1/e and 1/e 2 singularities between the I [2 — ► 2] and I [2 — > 3], 
including the counterterm / 2 ^ 2 [CT] defining the MS NLO parton distributions, is 
illustrated in the Table ^ below. 

The simple result of the whole procedure then is that the full cross section can be 
written as 

/ = I [2 -> 2 net] + I[2 -> 3 net] (15) 



7 



Term 


l/e 2 




finite 


/17b 3 [soft] 

A 2 ? 3 [soft] 
Ij^z [collinear] 
Ii 2* 3 [collinear] 
J 2 ^ 3 [finite] 
J 2 ^ 2 [Born] 

J 2 ^ 2 [HO] 
J 2 ^ 2 [CT] 


OA,B 
dl,2 

Ol,2; «A,B 


&A,B) C A,B 

1,(1). r 
°1,2) c l,2 

».( 2 ) J 

°a,b! c a,b; «a,b 

1,(2) 

& 1,2! Cl,2 

1,(1) .1,(2) . fe(l).fe(2) 
°A,B) °A,B) °1,2) °1,2 

^A,B 


X 
X 
X 
X 
X 
X 
X 



Table 1 : The singularity content of the different terms contributing to the cross section (^|) . 
The top panel specifies the type of the singularity and the symbols cij, b\ , b[ , Cj, di (i = 
A,B,1,2) represent the coefficient functions of the singularities of that type. Singularities 
between the different terms cancel whenever two identical symbols appear. For instance, the 
l/e 2 singularity in the term J 2 ^ 3 [soft] is cancelled by that contained in the NLO 2^2 
term 7 2 ^ 2 [HO]. The counterterm i~ 2 ^ 2 [CT] defines the NLO parton distributions in the MS 
scheme. Appearance of a finite contribution to the cross section (Q) in each term is also 
indicated by the crosses in the last column. 



where the 'net' contributions are obtained by removing all of the l/e and l/e 2 terms 
and letting e — > 0. Due to the infrared safety of this NLO algorithm the divergences 
originating from the higher order loop diagrams to 2 — > 2 scatterings and soft particle 



emission of 2 — > 3 scatterings cancel each other exactly |3(J |34 



2.2 Measurement function for the minijet Et 

We now turn to the detailed description of the implementation of the algorithm dis- 
cussed in the previous section. The main points have been addressed already in our 



previous work [3^] but for completeness, we wish to do this again in here. The quantity 
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we are first after, is the total Et carried by minijets into a chosen rapidity acceptance 
window in an average inelastic pp collision. Furthermore, in order for the perturbative 
treatment to be valid, the partonic collisions must be required to be sufficiently hard. 
All this has to be included in the measurement function which must be infrared safe 
by construction. 

For the computation of observable jets, the starting point would be to define the 
jet cone, i.e. the definition of when two nearly collinear final state partons are to be 
considered as one jet and when as as two separate jets. Such a measurement function 
is not applicable here but an acceptance region of a new type must be defined. 

In a hard scattering of partons in the NLO, we may have one, two, three or zero 
minijets in our rapidity acceptance region AY, defined in the (y, 0)-plane as 



AY = {{y,4>) : Vmin <V <Vr, 



< 4> < 2tt}. 



(16) 



In our case, AY will always be centered at y — 0, see Fig. [I]. As only massless partons 
are considered here, the transverse energy entering AY can be defined as a sum of the 
absolute values pri of the transverse momenta of those partons whose rapidities are 
within AY: 



E T = e(yi)p T i + e(y 2 )p T 2 + e(y 3 )pT3, 
where the step function e(?/j) is defined as in fl3fl , 

| 1 if yi G AY 
I otherwise. 



(17) 



2tt 




-0.5 



0.5 



Figure 1: The acceptance region Al" in the plane of rapidity y and azimuthal angle 
window Al" is centered around y = here. 



The 



Since we cannot solve full QCD, but have to resort to perturbative treatment, we 
must carefully specify which collisions are to be considered perturbative, i.e. hard 
enough. We define the perturbative collisions to be those with a large enough amount 
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(^> Aqcd) of transverse momentum produced, regardless of where the partons go in 
rapidity. For the different final state kinematics, this implies 



2 : 

3 : 



Pti + PT2 > 2p 
Pti + PT2 + PT3 > 2p , Po > A< 



-QCD 



(19) 



where the parameter po, which restricts the computation below, is a fixed external 
parameter which does not depend on AY. For 2 — > 2 processes, appearing both in 
LO and in NLO, where the scattered partons are back-to-back in the transverse plane, 
the condition above means that individual minijets, inside or outside Ay always have 
Pt > Po- F° r the 2 — > 3 processes, which appear only in NLO, a single parton may 
have also less than p transverse momentum as long as the total transverse momentum 
released in the partonic process exceeds 2p . This feature, which is also the main 



difference to the formulation of the same problem in Ref. [3~3|, will have interesting 
implications later on. We will investigate the dependence of the results on p Q and 
whether the NLO results remain stable relative to LO as po is taken into the few-GeV 
realm. To estimate the transverse energy produced in AA collisions at RHIC and LHC, 
additional phenomenology (see e.g ||22|| ) needs to be introduced. This phenomenology 
will be separately discussed in Sec. 4. 



The infrared safe measurement functions S 2 and S3 [30] can now be written down 
by combining the definitions of perturbativeness of the collisions and the definition of 
Et in AY. For 2 — > 2 scatterings, we define 

StipuPz) = ©(Pri +PT2 > 2p )5(E T - [e{y 1 )p T i + t(y 2 )pT2]) (20) 
and for 2^3 scatterings correspondingly 

S3(Pi,P2,P3) = ®(pti+Pt2+Pt3 > 2p )5(E T - [e(yi)p T i + t{y 2 )p T 2 + ^{yz)Prz])- (21) 

where O is a step function. 

It is easy to verify that these measurement functions indeed fulfill the requirements 
(0) and ([D for infrared safety. This is also intuitively clear, as the total Et is not 
sensitive to whether it is carried by one parton or two collinear partons, or by two 
partons of which one is at the soft limit. 

The measurement functions above define the semi-inclusive Et distribution of mini- 



jets in AF in pp collisions, introduced in LO in jr3[. Extending this now to NLO 
gives 



da 



dEn 



da 



VoAY 



dE- 



2^3 

PoAY 



2 ^ 2 da 
p ,AY dE T 

da 2 -* 2 r 

d[PSh d[PSh S2{fuP " ) + i d[PS] 



da 2 ^ 3 

3 irr ,m S 3 (Pl,P2,P3] 



'd[PS]i 



(22) 
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with the notation as in Eq. (|2|). 

The minijet Et distribution above is normalized to the integrated cross section 
(T iet(y/s,Po), as explained in [O. As a side remark, we note that although we do 
not explicitly consider the computation of a- ]Ct here, it is an infrared safe quantity to 
compute, i.e. the extension of its definition to NLO can be done. The main interest 
now, as already emphasized, is to study the behaviour of the first moment of the Et 
distribution (|22|) as a function of the parameter p which defines the extent of our 
perturbative calculation. Integrating the delta functions away in Eq. ( f22"D we get 



&qcv>\Et)ay,po = / dE T E T 

J 



do 



dEn 



PoAY 



°(E T y^ po + o{E T y^ 



where 



a{E T ) 2 - 2 



AY, P0 



a{E T )^ 



>d[PS] 3 

The measurement functions for the first EV-moment above are denoted by 



&(Pl,P2) 



PT2<$){PT2 > Po) 



(23) 

(24) 
(25) 

(26) 



e{yi)PTi + t{y 2 )p T 2 + e(y 3 )PT3 



©(Pti +PT2 +PT3 > 2p ), (27) 



where pri = PT2 in the former term and pri = |pt2 + Pt3| in the latter one. Naturally, 
also 5*2 and 5*3 fulfill the criteria (^) and (^[), which ensures that ct(Et)ay,p is a well- 
defined infrared safe quantity to compute. From the general form of the equations 
above, we notice that, by replacing the measurement functions S2 and S3 in Eq. |2| 
by and S3 defined above, we can exactly follow the formulation and solution of the 
problem as given in ref. [^1 an d summarized in the previous section. 

Finally, apart from the measurement functions themselves but related to the problem 
being infrared safely defined, the renormalization scale fiR in the strong coupling a s (fin) 
and the factorization scale fip m the parton distributions fi(x,fip) have to be chosen 
in such a way that the scales for the 2 — > 3 terms reduce to those for the terms with 
2 — > 2 kinematics in the soft and collinear limits. 



3 Numerical Evaluation 

The formulation of the minijet Et production thus exists. To study its implications 
and to estimate its usefulness a numerical study has to be carried out next. 
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Due to the simple structure of the two-particle phase space, one can analytically 
describe the kinematically allowed region with the cuts imposed by AY. The part 
consisting of two-particle final states is therefore easily evaluated with a high precision 
by using the available integration routines of the fortran NAG-library |}fj . 



The seven-dimensional phase space integral over the three-particle phase space re- 
duces to a six-dimensional one after the symmetry of the measurement function S3 
under the azimuthal angle <p2 is taken into account. In principle, these six-dimensional 
integrals could be treated similarly as the three-dimensional ones over the two-particle 
phase space, but due to the additional kinematical cuts introduced by the subtraction 
terms (the additional step functions in Eq. flT2])) this does not work out in practice. 
Rather one is forced to use fairly crude and quite time-consuming Monte Carlo integra- 
tion. We have used another NAG subroutine to take care of these terms. Fortunately, 
the major part of the cross sections comes from the Born level terms plus the finite 
parts of the subtraction terms having 2 — > 2 kinematics. The finite terms with 2^3 
kinematics are smaller and therefore their contribution to the total relative numerical 
error of the results is estimated to be rather small, below four per cent in any case. 

In the numerical evaluation of the integrands with 2 — > 3 we are able to use certain 
subroutines of the jet-program of Ellis, Kunszt and Soper J37|. Also the partonic 
book-keeping is taken care of by adopting the method and the permutation tables of 
different subprocesses used in the program of EKS. The number of flavors is fixed to be 
Nf = 4, and calculations were performed using two different sets of parton distribution 
functions (PDF): GRV94 ^ and a more recent set CTEQ5 |9f. For A QCD we take 
the value as quoted in the PDF set used. 

Before the numerical evaluation, also the renormalization and factorization scales 
have to be fixed. In principle, the dependence of the results on fiR and fip should be 
studied independently. In practise, however, the Monte Carlo integrals being tedious 
to solve numerically, we follow the common practise and choose these to be equal, 
/Ar = \i>F = A*, and study the dependence of the results on /i only. Regarding the 
scale itself, the most physical choice should reflect the perturbativeness of the collision, 
therefore we set [i proportional to the total transverse momentum produced in the hard 
process, regardless whether the partons are in AY or not, 

2^2: fi = iV M x -( m + p T2 ) = x p T (28) 

2^3: fi = N^x ~(p T1 + p T2 + Pts), 

where is a constant of the order of unity. This choice is infrared safe in the sense 
of Eqs. (0)-©, as is required for the exact cancellation of the divergent terms. We 
also note that with this scale choice one is able to include the new kinematical region, 
where two partons fall outside AY and one with small momentum (pt < Po) falls 
inside, contributing thus to the region Ej- < p not present in the computation of the 
perturbative 2^2 processes. 
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As there is some freedom in choosing the scale, other possibilities can also be thought 
of. For example in Ref. |$3| the total Et entering AY was considered. This, we think, 
leads to an ambiguity already at the 2 — > 2 level, as almost identical hard subprocesses 
are given very different scales depending whether both of the partons fly into AY or 
not. This issue is also clearly connected with the definition of the perturbativity of 
the partonic processes, i.e. which collisions are to be included in the computation and 
which excluded. For high-p r jets, of course, /i ~ pj? is the most natural scale choice, 
but in our case we would like to emphasize the semi-inclusive nature of the problem, 
and also that the minijet Et distribution is, unfortunately, not a direct observable, 
especially not in AA collisions. 



4 Results 

4.1 Dependence on po and on the PDF set 

The main goal, for the reason of the importance of the semihard region in initial Et 
production in AA collisions, is now to study the behaviour of ct(Et) as the parameter p Q 
controlling the perturbativity is taken from large scales into the few-GeV domain. On 



the basis of high-p^ inclusive jet computations [40] one may well expect fairly 

large NLO contributions: at scales of few tens of GeV, say, NLO improved results 
may well be larger by a factor of two as compared to LO ones. This largeness in the 
absolute magnitude of the correction terms does not imply that the perturbation series 
would be in a stall well out of the region of convergence, but rather is a signal of the 
possibility of some new element becoming manifest only at such high order. For the 
case of observable jets, the new element is the dependence of the NLO cross sections on 
the jet cone R. Another well-known example of precisely such a case is the Drell-Yan 



dilepton production where the NLO corrections bring in QCD as a new element [41 
In this sense, the question of the convergence of the perturbation series is moved to 
the behaviour of the next-to-next-to-leading order terms. 

In our case the new element the NLO definition of the problem introduces, is the 
appearance of a new kinematical region: for 2 — > 2 processes, the Et distribution is 
empty at scales E T < Po but not anymore so in the 2^3 case. The crucial question 
therefore concerns the behaviour of the results as po is decreased from, say, 10 GeV 
down to to 1... 2 GeV. If the NLO result is of a similar magnitude as at 10 GeV, as will 
be shown to be the case, the conclusion is that the NLO result is stable and pQCD 
can be appealed to at such small scales. Another interesting point, already studied in 



32], is the dependence of <j{Et) on the choice for the renormalization and factorization 
scale, and whether this behaviour substantially deviates from that in LO. 
To study the stability in the sense described above we define two ^-factors: 

K' = ^ K = ^ (29) 
LO" LO ' 1 ] 
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where LO' is the first term of the NLO expansion, i.e. it is the Born term evaluated 
with 2-loop a s and NLO parton distribution functions. Unprimed LO stands for lead- 
ing order result evaluated with 1-loop a s and LO parton distribution functions, this 
being the "truly" leading order result. These two i^-factors tell different tales: The 
unprimed K is the one commonly introduced in the LO computations to bring the 
elements of NLO in. The K' on the other hand measures roughly the difference be- 
tween two subsequent terms in the perturbation series, and possibly carries somewhat 
more information about the actual convergence of the perturbation series, since in K' 
the same parton distributions appear both in the nominator and denominator. Which 
one is more relevant, is somewhat of a matter of taste, though, and this is why both 
i^-factors are discussed. 



\y\ <0.5 ; GRV94 \y\ <0.5 ; CTEQ5 




Po [GeV] Po [GeV] 

Figure 2: (a) The first moment (j{Et) of the minijet Et distribution as a function of the 
minimum transverse momentum scale pq. Three curves for both RHIC, \fs = 200 GeV, and 
LHC, y/s = 5500 GeV, are shown: the solid lines are the full NLO results computed with the 
GRV94-HO parton distributions, the dashed lines are the LO results evaluated with 2-loop 
a s and NLO parton distributions, and the dotted lines in between are the truly LO results 
evaluated with 1-loop a s and GRV94-LO parton distributions. The scale choice for all points 
in this figure is = 1 in Eq. (^). (b) The same but computed with the CTEQ5M and 
CTEQ5L parton distributions. Labelling of the curves is same as in panel a. 



In Fig. |]a the results for a(E T ) computed with GRV94 parton distribution functions 
are shown as a function of the parameter p Q . The solid curves are the NLO results, 
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the dotted and dashed ones show the LO results computed in the two different ways 
described above. In Fig. we repeat the calculation by using the CTEQ5 parton 
distribution functions. The corresponding f^-factors (K and K') are collected in Table 
§ for some values of po from both figures. The immediate observation from the figures 
and the table is that the full NLO results deviate very systematically from the LO ones: 
Both K and K' increase as po —>■ 1 GeV, but the factor K appears to increase less than 
the factor K' . The steady increase in K' can be understood as a signal of the nearby 
borderline of pQCD. In a more practical sense the differences between LO and LO' arise 
at RHIC energies from the differences in a s and Aqcd arid at LHC energies, in addition 
to the previous, from the difference between the LO and NLO parton distributions. 

Based on the more stable behaviour of K, it may well be argued that if one trusts 
the perturbative computation at 10 GeV scale, one can equally well count on it also at 
2 or 3 GeV scales. Another observation is that the .ff -factors do depend on the PDF 
sets chosen. 



Po 


^RHIC 




TfiGRV(CTEQ5) 
LHC 


TfGRV(CTEQ5) 
LHC 


1 


4.8 (5.9) 


2.2 (3.4) 






2 


3.6 (4.2) 


2.2 (3.3) 


4.2 (4.9) 


1.5 (2.1) 


4 


3.0 (3.2) 


2.1 (3.0) 


3.4 (3.9) 


1.5 (2.2) 


6 


2.7 (3.0) 


2.0 (2.9) 


3.2 (3.5) 


1.6 (2.2) 


8 


2.4 (2.6) 


1.9 (2.7) 


3.0 (3.3) 


1.6 (2.2) 


10 


2.3 (2.5) 


1.8 (2.7) 


2.9 (3.1) 


1.6 (2.2) 



Table 2: Values of the two K- factors defined in the text for different values of pq. The 
numbers are shown for both RHIC and LHC energies and for two different PDF-sets, the 
GRV94 and CTEQ5. 

It should be emphasized once more, that here p is to be strictly regarded as an 
external parameter, on which the dependence of the results is studied. To form the 
initial in AA at RHIC and at the LHC, one has to decide upon at which value 
of po the ct{Et) should be taken. The expectation is that po ~ 1 ... 2 GeV, possibly 



depending on y/s and A |22| , but precise estimation of this value is beyond the pertur- 
bative approach in the sense that it requires introduction of further phenomenology or 
knowledge of non-perturbative physics involved. Regardless of the actual mechanism 
to determine the suitable value of po the basic question of applicability of pQCD in 
the few-GeV domain remains, and it is precisely this question that is currently being 
answered. 
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4.2 Dependence on the scale choice 



The dependence of <j{Et) on the renormalization/factorization scale /i, chosen accord- 
ing to in Eq. fl29|) , has essentially been presented in Ref. |32| but for completeness let 
us discuss this a bit further here. In Fig. H, we plot cr(E T ) at a fixed value p = 2 GeV 



as a function of the AL in Eq. 



The results in Fig. |2|(a) for this po are found at 



N/j, — 1 in the figure. 

At RHIC the results depend somewhat on the chosen scale. This behaviour is as 
what could be expected on the basis of the studies of cross sections of observable jets 
p9| , |40| , where it is more difficult to find a plateau as a function of fi/px towards 
smaller jet momenta. At the LHC <j(Et) exhibits a clear plateau around A" M ~ 1. 
Similar behaviour is observed already in both leading order curves, due to cancellation 
between the rise of gluon distributions and the decrease of a s with an increasing scale. 
Therefore the plateau cannot be attributed to a fast convergence of the perturbation 
series but the NLO results depend equally strongly (RHIC) or weakly (LHC) on the 
scale choice as the LO results do. Another consequence of this is that the A'-factor is 
close to constant near ~ 1. 
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Figure 3: The first moment ct(Et) as a function of = fJ-/J2i(PTi) (see Eq. p9|)) at a 
fixed value of po = 2 GeV. The labelling of the curves is same as in Figs. H and the GRV94 
parton distributions have been used. 
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4.3 A new kinematical region 



As noted already in our earlier work [04 , the NLO computation allows for a new 
kinematical region to be included. The existence of such a new contribution arises in 
the following way: in LO and in the loop-corrected NLO terms, where only 2^2 
processes are considered, the two final state partons are produced back-to-back with 
equal transverse momenta. The smallest amount of E? entering the acceptance region 
AF is therefore equal to po, and consequently the minijet ^-distribution is completely 
empty below Et = Po- In NLO, however, the 2^3 scatterings are included and the 
situation where one of the three partons has transverse momentum smaller than p but 
the other two have larger momenta, as required by the condition p T1 + p T2 + Pt3 > 2po 
and momentum conservation, may arise. Now, the two partons having transverse 
momenta larger than p may fall outside the acceptance window and only the third 
one with p? equal to a fraction of po enters AY , thus adding a contribution from the 
domain Ej- < p Q into the ^-distribution. 

To identify the contribution from the new domain at Et < Po from the results shown 
in Figs. we compute a{Ej) by performing the integration in Eq. (|23|) from Et = Po 
to t/s instead of integrating all the way from 0. To the measurement function ( |2"7j ) this 
adds a term Q(E T > p ). This is actually close to the approach of Ref. |33| (but still 
not the same, as our p is a constant independent of AY), where the new kinematic 
region is not included. The obtained results are shown in Fig 0, and the value of the 



resulting i^'-factor is comparable with the ones obtained in |33 . 

Comparing this figure with Fig. 0b, one observes that the i^-factors come down 
by a factor of two, hence enhancing the indications of the validity of the perturbative 
treatment. At the LHC energy, and at RHIC for po ~ 1 GeV, the NLO result shows 
practically no deviation from the truly LO one. Also, by comparison of the two figures, 
the slight growth of the K factors, especially relative to LO', can now be attributed 
to the inclusion of the new kinematical region at Et < Po- The conclusion thus is 
that a new kinematical domain at Et < Po not only exists, but it is also an important 
new feature in the extension of the problem to NLO. This also hints that if a conver- 
gence of the perturbation series is looked for in the semihard region, at least a NNLO 
computation should be studied. 

It is also interesting to note that inclusion of the new kinematical domain brings the 



perturbative approach towards the studies of the classical field configurations ^3 
[14] , [45|1 , where the dominant contribution is obtained from the process of two valence 
quarks exchanging a gluon and emitting a gluon. All such three-particle processes, 
and most importantly the purely gluonic ones, are now rigorously included in our 
computation, supplemented of course with the criterion for perturbativity, i.e. that in 
our approach the partons which scatter forward and backward outside AY, are required 
to carry large enough pt- 
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Figure 4: (a) The quantity u{Et) as a function of po, computed as in Fig. ||(b) but with 
an additional cut 6{Et > Eq = po))- Labelling of the curves is the same as in Figs. and 
CTEQ5 partons distributions are used. Note that the LO curves (dotted and dashed) are 
identical with those in Fig. ^(b). 

4.4 Dependence on AY 

So far all the results presented have been evaluated for the central rapidity unit, and 
it would certainly be desirable to know how the results change as the rapidity interval 
AF (centered around y — 0) is varied. This is shown in Fig. [5], where a(Ej) is plotted 
as a function of AY, keeping p at a fixed value of 2 GeV. The results are seen to be 
nearly linear in AY around AY = 1. 

In fact this is an important result as it increases confidence in computing the initial 



conditions from pQCD [22j in the following way: The average initial energy densities 

a! 

(30) 



in the central spatial region, z ~ 0, can be estimated through the Bjorken formula 

E$ A (b = 0,^, Po ,AY) 



nR 2 A T AY 

where the initial time r = 1/po and the spatial rapidity interval in the longitudinal 
volume element in the denominator has been taken to be AY . As seen in Fig. |, the 
initial transverse energy E^ A = Taa&{Et) behaves in LO essentially linearly in AY 
near the central rapidities and the AY dependence is thus cancelled out in the initial 
energy density of Eq. fl3"0"|). Based on the NLO computation of the cross sections of 
observable jets p9fl , it could be expected that AY plays the role of the jet cone R and 
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y&AY p Q =2.0 GeV ; CTEQ5 




0.2 1 2 

AY 

Figure 5: The quantity <j(Et) as a function of the rapidity acceptance window Ay with 
a fixed value of po = 2 GeV. The three lower curves are for RHIC and the three topmost 
curves for the LHC. As before, the solid curves stand for the NLO results, and the dotted 
and dashed ones for the LO results. 



the NLO results would therefore exhibit a different behaviour (~ R 2 at large R) in 
AF than the LO results. From the point of view of the local energy densities obtained 
through Eq. fl30|) this would be disastrous as would badly depend on the Ay chosen. 
As clearly seen in the figure, this is now not the case but the NLO results behave nicely 
in a similar way as the LO results do. This implies that the initial energy densities are 
not sensitive to the detailed size of the rapidity interval. One should also pay attention 
to the difference between the calculation of Et of ordinary jets and the calculation of 
Et in Ay from minijets considered in here: the jet calculation introduces an additional 
dependence on a resolution parameter, the jet cone radius R, which is simply absent 
in our considerations. 



4.5 Nuclear effects 

It is a well known fact that the distribution of partons in a free proton are different from 
those of a proton bound to a nucleus, fi/A(x, fi 2 ) = Rf(x,fj 2 )fi(x,fj 2 ) with Rf ^ 1. 
Regarding the behaviour of ct(Et), the last point of study here is the effect of an 
inclusion of nuclear effects Rf(x,Q 2 ) in the parton distributions. This is done by 
using the EKS98 parametrization [F2tJ which is based on a DGLAP analysis of the 
nuclear parton distributions, with constraints from the existing data from deep inelastic 



19 



I A scattering and Drell-Yan production in pA and also from momentum and baryon 
number conservation. There are definite sources of uncertainties we implicitly accept 
when using EKS98: First, the nuclear modifications have been shown to be independent 
of the PDF set chosen within a level of a few per cent error |2(J . Second, the EKS98 is 
based on a LO DGLAP evolution equations only but we apply it to the NLO PDFs as 
well. Finally, the nuclear gluon distributions themselves are still largely unknown at 
x<0.01, a region which is important for the few-GeV minijets at the LHC For RHIC, 
the situation is better since for 0.01<x<0.2 some constraints can be obtained from the 
Q 2 dependence of the ratio of the structure functions F 2 Sn /F 2 c measured by the New 
Muon collaboration in deeply inelastic lepton- nucleus scatterings pjf] . 
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Figure 6: (a) The a{Ex) as in fig, ^(a) but also nuclear shadowing has been taken into 
account via EKS98 parametrization PQ| . The computation has been carried our for a nucleus 
A=208 for both energies, (b) Comparison of the ct(Et) with nuclear shadowing (solid and 
dotted lines) and without shadowing (dot-dashed and dashed lines). The curves are shown 
for RHIC (lower four curves) and LHC (upper four curves). Of the four curves the upper 
pair is the full NLO result and the lower pair is the truly LO result. 



The validity region of the EKS98 does not quite extend down to \i — 1 GeV, so we 
limit the study of nuclear effects in the NLO study of a(E T ) to p between 2 and 10 
GeV. The results can be found in Fig. |6] below. In the first panel, the nuclear effects 
are included in all curves, the labelling of which is the same as in previous figures. The 
corresponding i^-factors are collected in Table [| Comparison of the numbers with 
those in parentheses in Table 2 indicates that the effects of the nuclear modifications to 
the -fT-factors can be neglected at RHIC, whereas for the LHC there is a slight increase 
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towards the region po ~ 2 GeV. The reason for this can be understood based on the 
second panel, where a comparison between the results with and without the nuclear 
effects is shown. For RHIC near po ~ 2 GeV, one is in the region where the nuclear 
effects of gluons stay very small as the cross-over region between shadowing (Rg < 1) 
and ant i- shadowing (Rg < 1) is probed (For 2^2, typically x ~ 2p /y/s). Towards 
Po ~ 1 GeV one would move into the shadowing region but the nuclear effects remain 
quite small in any case. At the LHC, the shadowing effects are larger, about 30% for 
the LO and 20% for the NLO results. The difference is due to the new kinematical 
region included: in the region Ex < Po, where typically y\ > 0.5 and yi < 0.5 and p^ is 
small, the values of x± and 22 are slightly larger than the typical in LO. The net 
loss in the gluonic luminosity due to shadowing is then altered less than in LO. 



Po 


K' 

^RHIC 


Krhic 


K' 


Klhc 
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4.2 


3.3 


5.5 


2.5 
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3.2 


3.0 


4.2 
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3.0 
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2.7 


2.7 


3.3 


2.2 


10 


2.6 


2.7 


3.1 


2.2 



Table 3: Values of the two iVfactors defined in Eq. (|29|) for different values of pq. The parton 
distribution functions are those of the CTEQ5-set and the nuclear shadowing is included 
through the EKS98 parametrization. 



5 From p to p sa J 

As described in the introduction, the goal of the present paper is to study the main 
uncertainties in the collinearly factorized minijet formalism to the initial production 
of transverse energy in ultrarelativistic nuclear collisions. Now that we have control 
over the NLO corrections in the quantity a [Ex) and have also verified that they do not 
present significantly divergent results at scales of few GeV, we may, as an application, 
turn to a determination of the actual value of the scale p$. As emphasized already 
in the introduction, this requires additional phenomenology. Here we will work in the 



framework of Ref. |22 



At certain scale po = p sat (^,A) the system of gluons produced into the central 
unit of rapidity will reach a density such that the effective total area of these gluons 
NAA(\/s,Psa.t) x Tr/Psat wn l be equal to the transverse area ttR\ in central collisions. 
Performing the computation at a scale po = p sat will effectively estimate the contribu- 
tions of all scales in the sense that gluons with px <C p sa t would contribute negligibly 
to the total Et- Hence, if p sat falls into the perturbatively computable region, above ~ 
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1 GeV, the initial energy densities can be obtained from the perturbative computation 
alone. For nuclei with A ~ 200 and collision energies of the order of 100 GeV this has 
been shown to be the case |2~2~ . 



To demonstrate the effects of the NLO terms and shadowing corrections, a determi- 
nation of the scale p sat is illustrated in the upper panel of Fig. [7], showing the number 
of minijets in a central AA collision of A = 208, Naa(P = 0,Po, y/s, AY) = T^(0) x 
2<7j et (po; y/s, AY) [|7| in the rapidity acceptance region AY. There is, however, a the- 
oretical caveat in computing the number of minijets, or rather 2aj e t(po> AY) = a(n), 
the first moment of the number distribution of minijets in AY to NLO. A strictly 
infrared safe definition of a measurement function for a(n) requires an introduction of 
an additional resolution scale, such as a screening mass, which defines when two nearly 
collinear partons within AY are to be counted as one and when as two. A detailed for- 
mulation of this is beyond the scope of the present paper, so let us simply simulate the 
NLO terms by taking the J\~-factors computed for a{Ej) and multiply the LO results 
for Naa by them. The value of po can then be determined by the intersection points of 
Naa(po) with the saturation curve Naa(Po) — Po^-a- ^ n the lower panel we show the 
results for the initial transverse energy E AA = Taa^(Et) with and without shadowing. 
The initial Et can then be read off from the different curves at each saturation scale 

PO = Psat- 

It is interesting to note, as pointed out already in |2^| , that at the saturation scale 
Po = Psat the computed E aa /Naa is very close to E/N = e^/r^ = 2.77] of an ideal gas of 
massless bosons, where the temperature T is obtained from e« = aT 4 with a = 167r 2 /30. 
From the point of view of energy per particle, a very rapid thermalization seems thus 
plausible. Recent studies [11], [[2], f||] indicate that this may indeed be the case, although 
early thermalization is still a subject of debate [j49[| . In any case, it is worth emphasizing 
that in our approach saturation and apparent thermalisation take place simultaneously. 

In fact, the early thermalisation could be used to get rid of the worry that the 
saturation approach involves a non-infrared safe quantity: one could first compute the 
thermalised initial parton number iVj(po) as a function of po from the infrared safe 
initial Et as E AA (p ) — > ei(po) — > Ti(p ) — > Sj(po) —> N^Po), and then use the obtained 
thermal iVj(po) i n the saturation criterion to determine p sat . 

If the system is thermal at p sat and if its further evolution is isentropic, the initial 
number of particles reflects the final multiplicity per unit rapidity quite closely, Nf w 
0.9NAA(j>sat, y/s) HU- Figure [7| can therefore be directly used to estimate the effects 
from the NLO contributions, parton distribution functions and shadowing to iVy If 
the NLO corrections are taken into account via the J\~-factor, then from the saturation 
criterion N A a = p'Lt^-A ^ follows, assuming that Naa ~ K/Psav that N A a ~ K 2 ^ n+2 \ 
Furthermore taking E T ~ Psat Naa one 

finds that E T ~ K 3 ^ n+2 \ From the Fig. 
we see that n = 3.0 . . . 3.3 at RHIC and n = 2.4 . . . 2.7 at LHC, depending whether 
shadowing is included or not. Let us take, say, K ~ 2 at RHIC and K ~ 1.5 at LHC. 
Then, using the smallest values for n quoted above, one finds that at RHIC the increase 
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Figure 7: Upper panel shows the determination of the scale p sa t- The saturation curve, 
Naa(po) = PqR\ i s shown by the dashed line labeled 'saturation'. 2 sets of four curves are 
shown. The upper four curves are for LHC energies and the lower ones are for RHIC. The 
solid lines are the LO results with GRV94 parton distributions multiplied by K (po) from the 
u(Et) analysis, and the dotted ones are the LO results with GRV94 without the ET-factor. 
The dashed lines stand for the LO results with CTEQ5 and the dot-dashed lines for the LO 
results with CTEQ5+EKS98, multiplied with respective K(p ) from the ct{Et) calculation 
with the same PDFs. The value of p sat can be read off in the intersection points of the 
'saturation'-curve with the other curves. The corresponding values of -Er(Psat) can be read 
off from the curves at the lower panel. 



with respect to the LO results in N and E T is 30% and 50% respectively. At LHC 
the increase is 20% iV and 40% in Et- This also demonstrates the importance of the 
inclusion of the NLO contributions 

It is more difficult to give the theoretical uncertainty in the saturation criterion it- 
self, due to the possible nonperturbative elements involved. The data coming now for 
the charged particle multiplicities from RHIC can in principle be used to constrain the 
remaining uncertainty in the saturation criterion. Once the scale p = p sat is nailed 
down for one collision energy, predictions of the final state Et within the hydrodynam- 
ical approach can be computed with the initial energy densities obtainable based on 
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the infrared safe NLO computation of E^ A . As mentioned in the introduction, a large 
decrease in Et is expected due to the work done by pressure. So far only longitudinally 
expanding QGP has been studied with minijet initial conditions [E3], the results from 
a transversally expanding hydrodynamic system will be reported soon |5C|. 



6 Conclusions 

We have presented in detail the theoretical basis, implementation and numerical results 
of a rigorous perturbative QCD calculation of the amount of transverse energy produced 
from minijets in AA or pp collisions into a rapidity acceptance region. The formulation 
is based on collinear factorization. In particular, the quantity <j{Et), the first moment 
of the minijet Et distribution, is computed and analysed. 

Theoretically, the key point is the formulation of an infrared safe measurement func- 
tion which gives the minijet Et distribution in AY and includes the definition of 
perturbative collisions. The singularities arising in the NLO terms are cancelled by 



applying the subtraction method []58, |3(|. The dependence of the results on the par- 
ton distribution functions, on the scale choice, on the width of the acceptance window 
and on nuclear effects in parton distributions have been studied in detail. Also, as an 
application, we have considered a way to fix the minimum transverse momentum scale 
Pq in the saturation model , and the uncertainties from this procedure to the initial 
conditions of the QGP in AA collisions at RHIC and LHC and also to the charged 
particle multiplicities in the final state. 

The analysis of the full NLO results leads us to conclude that the NLO corrections 
are large but well behaved in the sense that their relative magnitude as compared to 
the LO results is not diverging as subprocesses with smaller and smaller amount of 
transverse energy exchanged were taken into account. Furthermore, the contributions 
from two different regions Et > Po and from Et < Po were studied. The latter one 
contributes only in NLO, through configurations where one small-pr parton enters AY 
and the other two partons with larger transverse momenta fall outside the acceptance 
window. When the contribution from the new domain Et < Po is excluded, the NLO 
corrections are observed to reduce by roughly a factor of two, and the stability of the 
results against the LO to increase. If the new domain is included, the comparison 
of the full NLO results against the LO numbers is not strictly one-to-one, since the 
new kinematical region is not present in LO at all. A consistent comparison in this 
case can therefore only be performed between NNLO and NLO contributions, as then 
contributions from all values of pt would be present. Our conjecture is that in that 
case suppression of NNLO terms will be observed. 

Another important and beforehand not at all obvious result is that the initially 
produced Et in AA collisions depends, contrary to what one might expect based on 
the jet cone dependence of high-py jets, practically linearly on the chosen rapidity 
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window AY. This has the important consequence that the local energy density does 
not depend on AY, and a possible source of ambiguity in estimation of the initial 
energy densities for the QGP in high energy AA collisions is removed. 

Finally, in connection to other models, in particular to the McLerran-Venugopalan 
model for gluon production through classical fields [Q, the inclusion of the new kine- 
matical region in our computation should bring these two approaches closer to each 
other. 
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